#!/usr/bin/env python
# coding: utf-8

# In[1]:


import os
import tensorflow as tf
os.environ["CUDA_VISIBLE_DEVICES"] = '3' #use GPU with ID=0
config = tf.ConfigProto()
config.gpu_options.per_process_gpu_memory_fraction = 1 
config.gpu_options.allow_growth = True 
sess = tf.Session(config=config)


# In[2]:


###### import numpy as np
import pandas as pd
import keras
# import tensorflow as tf
from datetime import datetime
from sklearn import datasets, linear_model
from matplotlib import pyplot as plt
from matplotlib import cm
from keras.layers.core import Dense, Activation, Dropout
from keras.layers.recurrent import LSTM
from keras.layers import BatchNormalization, Input, Embedding, Concatenate, Conv1D, MaxPooling1D, Flatten, merge
from keras.layers import merge, Concatenate, Permute, RepeatVector, Reshape,Multiply
from keras.models import Sequential, Model
import keras.backend as K
import statsmodels.formula.api as smf
from keras.callbacks import ModelCheckpoint
from keras.preprocessing.text import Tokenizer
from keras.preprocessing.sequence import pad_sequences

# prevent tensorflow from allocating the entire GPU memory at once
config = tf.ConfigProto()
config.gpu_options.allow_growth=True
sess = tf.Session(config=config)


# In[3]:


import pandas as pd
import math
import keras
import numpy as np


# In[4]:


data=pd.read_csv('station_all.csv')
# print(data.head)


# In[5]:


station=data.iloc[:,11:]#天气+十天发电量
trend=data['trend']
std=data['std']
train=data['detrended']#y_train
# print(train)
len_train=math.ceil(len(station)*0.7)
len_val=math.ceil(len(station)*0.2)
len_test=math.ceil(len(station)*0.1)

#print(std.head)


# # 划分训练集，验证集，测试集

# In[25]:


#最后每个训练集的最后一个
i_train=len_train
i_val=len_val+len_train
i_test=len_test+len_val

#训练数据
station_train = station.iloc[:i_train,:]
y_train=train[:i_train]#y_train==detrend
trend_train = trend[:i_train]
std_train=std[:i_train]

print(np.shape(station_train))#(433991,14)
print(np.shape(y_train))#(433991,)

print(np.shape(station_val))
print(np.shape(y_val))

#验证数据
station_val = station.iloc[i_train:i_val]
y_val = train[i_train:i_val]
trend_val  = trend[i_train:i_val]
std_val =std[i_train:i_val]

#测试数据
station_test = station.iloc[i_val:]
y_test=train[i_val:]
trend_test = trend[i_val:]
std_test=std[i_val:]

# len_train = 365*2#2017-2019
# len_val = 365*3#
# len_test = -1 # 2020-现在


# # 评价函数

# In[7]:


def compute_error(trues, predicted):#    回归评价函数：MAE、RAE、RMSE、RRSE、MAPE、R2
    corr = np.corrcoef(predicted, trues)[0,1]#计算皮尔逊积矩相关系数 r
    
#     皮尔逊相关系数是一种度量两个变量间（线性）相关程度的方法。它是一个介于 1 和 -1 之间的值，
#     其中，1 表示变量完全正相关， 0 表示无关，-1 表示完全负相关。
   
    mae = np.mean(np.abs(predicted - trues))#平均绝对值误差（MAE）
#      """
#     numpy.mean(a, axis, dtype, out，keepdims )   求取均值
#     经常操作的参数为axis，以m * n矩阵举例：
#     axis 不设置值，对 m*n 个数求均值，返回一个实数
#     axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
#     axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵
#     操作对象为：数组或矩阵
#     数组：
#     a = np.array([[1, 2], [3, 4]])
#     >>>np.mean(a, axis=0) # axis=0，计算每一列的均值
#     array([ 2.,  3.])
#     >>> np.mean(a, axis=1) # 计算每一行的均值 
#     array([ 1.5,  3.5])
#     """
    rae = np.sum(np.abs(predicted - trues)) / np.sum(np.abs(trues - np.mean(trues)))
#      """
#     sum(a, axis=None, dtype=None, out=None, keepdims=np._NoValue)
#     a：用于进行加法运算的数组形式的元素
#     axis的取值有三种情况：1.None(在默认/缺省的情况下），2.整数， 3.整数元组。
    
#     如果axis取None，即将数组/矩阵中的元素全部加起来，得到一个和。
#     如果axis为整数，axis的取值不可大于数组/矩阵的维度，且axis的不同取值会产生不同的结果。
#     如果axis为整数元组（x，y），则是求出axis=x和axis=y情况下得到的和。
#     """     
    rmse = np.sqrt(np.mean((predicted - trues)**2))#均方根误差（RMSE）
    rrse = np.sqrt(np.sum((predicted - trues)**2) / np.sum((trues - np.mean(trues))**2)) #根相对平方误差（RRSE）
    mape = np.mean(np.abs((predicted - trues) / trues)) * 100#平均绝对百分误差（MAPE） 
    r2 = max(0, 1 - np.sum((predicted - trues)**2) / np.sum((trues - np.mean(trues))**2))#R-平方  ,max 返回给定参数的最大值
    return corr, mae, rae, rmse, rrse, mape, r2


def compute_error_filtered(trues, predicted, filt):#计算错误已过滤   filt过滤器
    trues = trues[filt]
    predicted = predicted[filt]
    corr = np.corrcoef(predicted, trues)[0,1]
    mae = np.mean(np.abs(predicted - trues))
    mse = np.mean((predicted - trues)**2)
    rae = np.sum(np.abs(predicted - trues)) / np.sum(np.abs(trues - np.mean(trues)))
    rmse = np.sqrt(np.mean((predicted - trues)**2))
    r2 = max(0, 1 - np.sum((trues-predicted)**2) / np.sum((trues - np.mean(trues))**2))
    return corr, mae, rae, rmse, rrse, mape, r2


# In[8]:


NUM_LAGS=10#滞后时间


# # FN+station+W

# In[27]:


def build_model(num_inputs, num_lags, num_preds):#  MLP  多层感知机
    #input_lags = Input(shape=(num_lags,))
    input_lags = Input(shape=(14,))
    x = input_lags

  

    x = BatchNormalization()(x)
#     print('输出的x:',x)
#BN算法可以省去这两项或者只需要小的L2正则化约束。原因，BN算法后，参数进行了归一化，
#     原本经过激活函数没有太大影响的神经元，分布变得明显，经过一个激活函数以后，
#     神经元会自动削弱或者去除一些神经元，就不用再对其进行dropout
#     另外就是L2正则化，由于每次训练都进行了归一化，就很少发生由于数据分布不同导致的参数变动过大，带来的参数不断增大。
#     可以吧训练数据集打乱，防止训练发生偏移
#     在卷积中，会出现每层卷积层中有（L）多个特征图。AxAxL特征矩阵。我们只需要以每个特征图为单元求取一对γ、β。
#     然后在对特征图进行神经元的归一化。

#定义基本的网络层
    x = Dense(units=100, activation="tanh", kernel_regularizer=keras.regularizers.l2(0.05))(x)
#     units:该层有几个神经元
#     activation:该层使用的激活函数
#     use_bias:是否添加偏置项
#     kernel_initializer:权重初始化方法
#     bias_initializer:偏置值初始化方法
#     kernel_regularizer:权重规范化函数
#     bias_regularizer:偏置值规范化方法
#     activity_regularizer:输出的规范化方法
#     kernel_constraint:权重变化限制函数
#     bias_constraint:偏置值变化限制函数
   
    x = Dropout(0.5)(x)#这防止过拟合，按照0.5的比例去掉一些神经元。
    preds = Dense(units=num_preds)(x)
    
    model = Model(input_lags, preds)
    model.compile(loss="mse", optimizer="adam")#loss被称为损失函数，optimizer为优化器（编译神经网络的两个参数）
    
    return model, input_lags, preds
checkpoint = keras.callbacks.ModelCheckpoint("weights2.best.hdf5", monitor='val_loss', verbose=0, save_best_only=True, mode='min')
# 该回调函数将在每个epoch后保存模型到filepath,filename：字符串，保存模型的路径
# monitor：需要监视的值，通常为：val_acc 或 val_loss 或 acc 或 loss
# verbose：信息展示模式，0或1。为1表示输出epoch模型保存信息，默认为0表示不输出该信息
# save_best_only：当设置为True时，将只保存在验证集上性能最好的模型
# mode：‘auto’，‘min’，‘max’之一，在save_best_only=True时决定性能最佳模型的评判准则，例如，当监测值为val_acc时，
#     模式应为max，当检测值为val_loss时，模式应为min。在auto模式下，评价准则由被监测值的名字自动推断。
# save_weights_only：若设置为True，则只保存模型权重，否则将保存整个模型（包括模型结构，配置信息等）
# period：CheckPoint之间的间隔的epoch数
model, input_lags, preds = build_model(1, NUM_LAGS, 1)
model.fit(
    #np.concatenate([lags_train], axis=1),#前两年数据
    station_train,
    y_train,
    batch_size=64,
    epochs=500,
    #validation_data=(np.concatenate([lags_val], axis=1), y_val),
    validation_data=(station_val, y_val),
    callbacks=[checkpoint],
    verbose=0)   
print ("共多少次 ", len(model.history.history["loss"]))
print ("最好的情况    ", np.argmin(model.history.history["loss"]), "   Best:", np.min(model.history.history["loss"]))
print ("验证集最好的情况", np.argmin(model.history.history["val_loss"]), "   Best:", np.min(model.history.history["val_loss"]))

model.load_weights("weights2.best.hdf5")

preds_lstm = model.predict(station_test)
# print(preds_lstm.shape,std_test.shape,trend_test.shape)
preds_lstm = preds_lstm[:,0] * std_test + trend_test
y_true = y_test * std_test + trend_test
corr, mae, rae, rmse, rrse, mape, r2 = compute_error(y_true, preds_lstm)
print ("MAE:  %.3f\tRMSE: %.3f\tR2:   %.3f" % (mae, rmse, r2))


# # LSTM+station+W

# In[ ]:

print('LSTM+station+W')


def build_model_weather(num_inputs, num_lags, num_feat, num_preds):
    input_lags = Input(shape=(num_lags, 1))
    input_weather = Input(shape=(num_feat,))

    x_lags = input_lags
    x_lags = LSTM(20,
                  kernel_regularizer=keras.regularizers.l2(0.1),
                  recurrent_regularizer=keras.regularizers.l2(0.1),
                  return_sequences=False)(x_lags)
    x_lags = BatchNormalization()(x_lags)

    feat = Concatenate(axis=1)([x_lags, input_weather])

    feat = BatchNormalization()(feat)
    # preds = Dense(units=num_preds)(feat)
    preds = Dense(units=num_preds, kernel_regularizer=keras.regularizers.l2(0.01))(feat)
    preds = Activation("linear")(preds)

    model = Model([input_lags, input_weather], preds)

    rmsp = keras.optimizers.RMSprop(lr=0.001, rho=0.9, epsilon=1e-08, decay=0.0)
    # model.compile(loss="mse", optimizer=rmsp)
    model.compile(loss="mse", optimizer="rmsprop")

    return model, input_lags, preds


checkpoint = ModelCheckpoint("weights2.best.hdf5", monitor='val_loss', verbose=0, save_best_only=True, mode='min')
station_test = station_test.values[:, 4:]
station_train = station_train.values[:, 4:]
weather_train = np.array(weather_train)
weather_test = np.array(weather_test)
station_val = station_val.values[:, 4:]
weather_val = np.array(weather_val)
model, input_lags, preds = build_model_weather(1, NUM_LAGS, 4, 1)
model.fit(
    # lags_train,
    # np.concatenate([lags_train, weather_feats_train[:,sel]], axis=1),
    [station_train[:, :, np.newaxis], weather_train[:, :4]],
    y_train,
    batch_size=64,
    epochs=1000,
    validation_data=([station_val[:, :, np.newaxis], weather_val[:, :4]], y_val),
    callbacks=[checkpoint],
    verbose=0)

print ("共多少次 ", len(model.history.history["loss"]))
print ("最好的情况    ", np.argmin(model.history.history["loss"]), "   Best:", np.min(model.history.history["loss"]))
print ("验证集最好的情况", np.argmin(model.history.history["val_loss"]), "   Best:", np.min(model.history.history["val_loss"]))

model.load_weights("weights2.best.hdf5")

preds_lstm = model.predict([station_test[:, :, np.newaxis], weather_test[:, :4]])
preds_lstm = preds_lstm[:, 0] * std_test + trend_test
y_true = station_test * std_test + trend_test
corr, mae, rae, rmse, rrse, mape, r2 = compute_error(y_true, preds_lstm)
print ("MAE:  %.3f\tRMSE: %.3f\tR2:   %.3f" % (mae, rmse, r2))
print (mape)


